Single-cell transcriptomic analysis of honeybee brains identifies vitellogenin as caste differentiation-related factor

Summary The honeybee (Apis mellifera) is a well-known eusocial insect. In honeybee colonies, thousands of sterile workers, including nurse and forager bees, perform various tasks within or outside the hive, respectively. The queen is the only fertile female and is responsible for reproduction. The queen and workers share similar genomes but occupy different caste statuses. We established single-cell transcriptomic atlases of brains from queens and worker subcastes and identified five major cell groups: Kenyon, optic lobe, olfactory projection, glial, and hemocyte cells. By dividing Kenyon and glial cells into multiple subtypes based on credible markers, we observed that vitellogenin (vg) was highly expressed in specific glial-cell subtypes in brains of queens. Knockdown of vg at the early larval stage significantly suppressed the development into adult queens. We demonstrate vg expression as a "molecular signature" for the queen caste and suggest involvement of vg in regulating caste differentiation.

Highlights scRNA-seq revealed distinct gene expression in the brains of queens and workers Vitellogenin (vg) may represent a "molecular signature" of the queen caste Knockdown of vg at early larval stage suppressed development into adult queens Vg may be involved in regulating caste differentiation

INTRODUCTION
The honeybee is a well-known eusocial insect. In a colony of honeybees, the queen is primarily responsible for reproduction, whereas the workers are sterile and perform tasks required for the colony's growth, maintenance, and defense (Johnson, 2010;Page et al., 2012). Workers usually perform tasks based on an agerelated (behavioral-developmental) process (Beshers and Fewell, 2001;Calderone, 1998;Robinson and Huang, 1998). In the first two to three weeks following eclosion, an adult worker is responsible for the nursing tasks (e.g., caring for the brood and queen); the worker subsequently switches to other specialized work for a few days, including building and guarding within the hive. Eventually, the role of the worker changes to foraging for foods such as nectar, pollen, and water, for the rest of its life (Ament et al., 2010;. This behavioral switch is common but not rigid. Instead, honeybees are sensitive to environmental changes and duties are adjusted, reflecting typical phenotypic plasticity (Duncanm et al., 2020;Maleszka, 2018;Sommer, 2020).
Queens and workers develop from diploid fertilized eggs, and therefore share similar genomes. However, they exhibit tremendous differences regarding behaviors and physiologies. The process involved in formation of the two castes is termed caste differentiation, which largely depends on feeding conditions during the early larval stage (within the first 96 h) (Kamakura, 2011;Wang and Li-Byarlay, 2015;Weaver, 1966). However, the molecular basis of the establishment and regulation of caste differentiation remains largely unknown. Previous studies proposed that nutritional conditions during early development profoundly influence the individual's epigenetic status (e.g., DNA methylation (Foret et al., 2012;Kucharski et al., 2008;Lyko et al., 2010), RNA m6A modification (Wang et al., 2021), and histone acetylation (Spannhoff et al., 2011)) and, thereby, determine the caste fate of the honeybee. In addition, nutrition-related signaling pathways, such as the insulin/insulin-like signaling (Nijhout and McKenna, 2018;Wang et al., 2013;Wheeler et al., 2006) and target of rapamycin nutrient-sensing pathways (Mutti et al., 2011;Patel et al., 2007;Wheeler et al., 2014), have been reported to participate in caste differentiation regulation.
Previous studies investigating neuronal cells in the bee brain suggested that there are remarkable differences in the transcriptomes (Behura and Whitfield, 2010;Liu et al., 2011Liu et al., , 2012Tholken et al., 2019), proteomes (Hernandez et al., 2012), and methylomes (Foret et  iScience Article queens. To better understand the difference between queens and workers, we established high-quality singlecell transcriptomic atlases of female honeybee brains. We identified a specific cell type in the brains of queensensheathing glial cells that highly express vitellogenin (vg). The expression of vg in a subset of glial cells serves as a "molecular signature" for the queen caste. Vitellogenin is a conserved yolk precursor protein exhibiting pleiotropic functions in honeybees. It has been reported that honeybee vg is associated with many central biological processes, including lifespan, reproduction, and immunity (Amdam et al., 2004;Brandt et al., 2005;Guidugli et al., 2005). To test whether vg contributes to caste differentiation, we performed RNAi (RNAi) assays and found that the knockdown of vg at the early larval stage significantly suppressed the development from larvae to the adult queen when reared under high-nutrition conditions. Therefore, in addition to its previously identified role in behavioral transitions between nurses and foragers (Ihle et al., 2010;Nelson et al., 2007), our findings suggest that vg is involved in regulating caste differentiation.

Single-cell transcriptomic atlases of honeybee brains
To identify caste-specific cell types among female honeybee brains, we employed 10X Genomics technology to conduct single-cell RNA-seq analysis. For this assay, we collected pools of four to five female honeybee brains from queens, foragers, and nurses, respectively. For each caste or subcaste, we sequenced two independent biological replicates, targeting over 10,000 cells per replicate. The 10X single-cell libraries were constructed, followed by high-throughput sequencing ( Figure 1A). Following the recommended standard procedures, FASTQ files were analyzed using the Cell Ranger pipeline to generate featurebarcode matrices. A total of 115,169 cells were obtained after data filtration using the Seurat R package with three parameters: the nCount (number of unique molecular identifiers), nFeature (number of genes), and percent.mt (portion of mitochondrial genes) ( Figure S1A). These cells included 40,186 (34.89%), 38,149 (33.12%), and 36,834 (31.98%) cells from the brains of foragers, nurses, and queens, respectively (Table 1). To validate the sensitivity and integrity of the single-cell dataset, we also performed bulk RNA-sequencing in parallel using corresponding samples. The comparison confirmed that the filtered single-cell data were highly correlated (r R 0.7) to our bulk sequencing data (Figures S1B and S1C).
In addition to the neurons, we found that cells expressing LOC410151 (rx2), the homolog of repo, were recognized as glial cells (Edwards and Meinertzhagen, 2010;Shah et al., 2018). As a result, 7,268 ($6.3%) glial cells were detected in our datasets, similar to the proportion found in Drosophila ($6.4%) (Davie et al., 2018). Moreover, we also observed a few cells ($0.3%, cluster 41) expressing the hemocyte markers LOC411597 (hml) and LOC551684 (Fer2LCH) that have also been found in the brains of fruit flies and ants (Davie et al., 2018;Sheng et al., 2020). Owing to the lack of marker genes enabling further subdivision, we termed the remaining 12 clusters "undefined neurons", because elav could be detected in these groups.
To compare differences in the cellular composition of brains derived from different castes, we counted the relative abundance of each cell category in the three female honeybee groups (Figures 1F, S2A and S2B). As shown in Figure 1F, the neuron compositions of queens and foragers exhibited similar patterns, and the ratios of KCs and OPNs were higher than those in nurses. Moreover, consistent results were obtained between the two biological replicates for each tested group.

Neuron classification using neurotransmitters and neuropeptides
During annotation, we observed that various genes associated with neurotransmitters and receptors were differentially expressed in the neuron subtypes. It has been demonstrated that neurotransmitters, including acetylcholine (ACh), glutamate, GABA, and monoamine molecules (e.g., octopamine and dopamine), are important in mediating distinct pathways involved in a variety of cognitive processes (Blenau and Baumann, 2001). On analyzing our datasets, we identified cholinergic (VAChT, ChAT), glutamatergic (Gad1, VGAT), and GABAergic (VGAT) neurons ( Figure 2A), based on established markers that participate in the synthesis or vesicle transport of neurotransmitters (Allen et al., 2020;Brunet Avalos et al., 2019), accounting for 40, 8, and 15% of the total neurons, respectively. Vmat has been used to label monoaminergic neurons (Allen et al., 2020), and we found that Vmat-positive cells ($10% of neurons) could be further divided into dopaminergic (ple), octopaminergic (Tdc2), serotonergic (DAT), and histaminergic (TDH) neurons (Figures 2A and 2B). Amongst the monoaminergic neurons, there were remarkably more dopaminergic neurons than the other three subtypes.
Previous studies reported the co-existence of more than one neurotransmitter in the same neuron in Drosophila adult and larval brain cell atlases (Brunet Avalos et al., 2019). To ascertain if this also applied to honeybee brains, we evaluated all possible combinations between the four markers VAChT, Gad1, VGlut, and Vmat. As shown in Figures 2C and 2D, nearly half of the neurons (51.9%) expressed one neurotransmitter, whereas the other neurons expressed two or more neurotransmitters. The statistical analysis revealed no significant differences in the neurotransmitter expression patterns among the three female groups ( Figure S3A). Notably, the frequency of dual (17.3%) or multiple (3.7%) neurotransmitter-expressing neurons in the honeybee brains was significantly higher than that in fruit flies (10 and 1.8%, respectively) (Brunet Avalos et al., 2019), revealing the complexity of the neural communication network in honeybee brains.
Neuropeptides represent another type of essential messenger molecule for communication between neurons. Although numerous putative neuropeptides were previously observed in honeybees, the role of these iScience Article peptides is poorly understood (Nassel and Homberg, 2006). We curated a panel of 27 reported honeybee neuropeptides and examined their expression at the single-cell level, and the heatmap in Figure 2E shows the cell-type-specific expression patterns of these neuropeptides in three female groups. Compared with results analyzed by bulk RNA-seq data ( Figure S3B), single-cell RNA-seq data showed a detailed differential expression in each cell type. For example, tachykinin and FMRFamide were highly expressed in mushroom bodies and optic lobes ( Figure S3C), which is consistent with previously reported findings (Schurmann and Erber, 1990;Takeuchi et al., 2004). Orcokinin, CTL9, and ILP-2 were enriched in Kenyon cells, glia, and hemocytes, respectively. By comparing the datasets from the three female groups, we found that although most peptides exhibited a similar expression pattern, a few displayed caste-specific patterns. For instance, Nus1 and apidaecin 1 were highly expressed in the hemocytes and glial cells of queens, respectively, whereas Crz expression was specifically present in worker subcastes ( Figure 2E).

Reclustering of Kenyon cells
Because KCs and glial cells are two major cell types in the honeybee brain, we then focused on elucidating their transcriptomic features. Kenyon cells of two major classes (class-I and class-II) located in the mushroom bodies are vital for sensory information processing, learning, and memory (Farris, 2005;Mobbs, 1982). Previous studies involving imaging analysis led to the class-I KCs being grouped into three subtypes: large (lKCs), middle (mKCs), and small (sKCs) (Suenami et al., 2018). By analyzing our single-cell seq datasets, we confirmed the existence of lKCs and sKCs with the specific markers mblk-1 and E74 (Suenami et al., 2016(Suenami et al., , 2018, respectively. Moreover, significant clustering of a KC population positively expressing FoxP (FoxP KCs) (Schatton and Scharff, 2017) was confirmed by our findings ( Figures 3A, and 3B). In addition, class-II KCs failed to be defined because of a lack of information on suitable marker genes for bees. It is noteworthy that the three female castes shared a repertoire of closely related KC subtypes, although the proportion of all KCs in nurses was less than those in foragers and queens ( Figure 3C).
We subsequently employed FindMarkers to analyze the differentially expressed genes (DEGs) in each KC subtype (Table S11) and found that the gene expression profiles of the lKCs, sKCs, and FoxP KCs between queen and worker castes were comparable on the volcano plot ( Figures 3D to 3F). The Venn plots in Figure S4A summarize the number of DEGs in the two groups for each cell subtype. It appeared that there were considerably fewer DEGs in each KC subtype between nurses and foragers than between foragers and queens, which was in agreement with the bulk RNA-seq results ( Figure S4B). However, the DEGs in these comparisons displayed no significant enrichment of gene ontology items.

Reclustering of glial cells
Given that glial cells are essential for the functioning of the nervous system, we performed glial cell reclustering analysis by employing the marker genes used for classifying different subtypes of glia in Drosophila (Allen et al., 2020;Parker and Auld, 2006). As shown in Figures 4A and 4B, glia reclustering revealed four subtypes: surface glia (Tret1, vkg), astrocytes (ebony, GlnS), ensheathing glia (Idgf4, Tsf1), and cortex glia (wrapper, zyd). By calculating the proportion of each cell type, we identified different compositions among the three female groups ( Figure 4C). In particular, the relative ratio of ensheathing glia was higher in queens' brains than in workers' brains. Ensheathing glial cells have been reported to function as phagocytes in the Drosophila brain (Doherty et al., 2009), and are also associated with the aging process (Sheng et al., 2020). We observed that nurse bees had relatively more ensheathing glia than foragers. Recent studies in Harpegnathos reported the expansion of ensheathing glia in the reproductive caste (Sheng et al., 2020). Hence, a high ratio of ensheathing glia could be a biomarker for the queen iScience Article caste. In addition, DEGs in each glial cell subtype were depicted in volcano plots ( Figures 4D-4G, Table S11), and the number of DEGs was summarized in Venn plots ( Figure S4C).

Ensheathing glia expressed vg in a queen-specific manner
Having observed the high number of ensheathing glia in the brains of queens, we sought to analyze the differential gene expression between the ensheathing glia of queens and workers. As shown in Figures 4G and S4C, 111 genes were more highly expressed in the ensheathing cells of queens compared to those of workers (nurse and forager bees). In particular, ensheathing cells in queens displayed a log2fold-change exceeding three in vg expression compared to the average expression levels of vg in worker ensheathing cells ( Figure 5A). To confirm this observation, we performed fluorescence in situ hybridization (FISH) assays, using multiple probes including vg, idgf4, and tsf1. In this assay, idgf4 and tsf1 were used to specifically mark the ensheathing glial cells in adult brains. As shown in Figures 5B-5D, vg transcripts were highly expressed in ensheathing glial cells in queens' brains compared to that in the ensheathing glial cells in workers' brains. Although the expression of vg varied in honeybees at different ages, we confirmed that vg was significantly more highly expressed in the queen caste rather than in an age-dependent manner by using qPCR analysis ( Figure S5A). As such, our results suggest that vg-positive ensheathing cells represent a specific subtype of glia in the queen brain and are likely to be a biomarker for the queen caste. Moreover, based on our findings, we then investigated whether vg contributes to caste differentiation.
Indeed, our bulk RNA-seq analysis also showed that vg was one of the most differentially expressed genes in the brains of queens and workers (Table S12). These results indicate that high expression of vg may play more important roles in queens compared to that in workers. Therefore, examining whether vg is involved in caste differentiation is an interesting issue that warrants exploration. Previous studies suggested that knocking down vg expression in adult workers affects the behavioral transition between nurses and foragers (Nelson et al., 2007). Given that caste differentiation is determined during the early larval stage, we knocked down vg throughout the larval stages by conducting RNAi on larvae kept in an in-dish rearing system ( Figure 5E). In this assay, we transferred the first instar larvae to a 24-well plate and fed them with adequate high-nutrition food compounds for development. For RNAi treatment, we packaged double-strand (ds)RNAs within a previously described synthetic nanomaterial (Jian et al., 2009;Lynn and Langer, 2000), which was added to the bee food ( Figure S5B). The bees were fed with dsRNA-treated food until the pupal stage. Pupae were then transferred to another 24-well plate ( Figure S5C), and the emerged adults were used for phenotypic analysis.
There are numerous differences in the body structure of queens and workers (Ament et al., 2010). Queen bees can easily be recognized by specific morphological features, particularly birth weight and body size. In our indish-rearing system, three groups of adults with different body phenotypes developed ( Figure 5F): (1) over 22% of adults appeared similar to natural queens (birth weight >190 mg), and we termed this group the "queen-like" bees; (2) $75% of adults were lighter than queens but heavier than workers (birth weight 120-190 mg), and we termed them "intermediate" bees; and (3) $3% of adults had body weights similar to those of workers (birth weight <120 mg), and we called this the "worker" phenotype. The worker phenotype was also confirmed at gene expression level using DEGs from bulk RNA-seq ( Figure S5D to S). As shown in Figure 5G, GFP dsRNA treatment did not affect the ratios of the aforementioned types of adults in our cultured system. However, vg knockdown significantly reduced the percentage of queen-like bees and increased the ratio of workers in tested populations. Tor is a known factor involved in caste differentiation (Mutti et al., 2011). Interestingly, we found that Tor knockdown displayed similar results as vg knockdown ( Figure 5G). Moreover, to test the role of vg at the early larval stage, we performed another RNAi assay in which vg was knocked down at the late larval stage (after 96 h) when the queen fate had already been established (Weaver, 1966). We did not observe any effects in this assay ( Figure S5T), indicating that vg contributes to caste differentiation during the early larval stage. Therefore, our results reveal that vg knockdown at the early larval stage can affect the development trajectory of queens, highlighting the importance of vg expression in bee development.  Figure S4 and iScience Article DISCUSSION Caste differentiation is a hallmark of a social insect community. The honeybee is an excellent model organism for studying the mechanisms underlying determination and regulation of the caste fate of queens. Previous studies proposed that nutritional conditions in coordination with epigenetic regulation and nutritionrelated signaling pathways control caste differentiation (Kucharski et al., 2008;Mutti et al., 2011;Spannhoff et al., 2011;Wang et al., 2021). Given that the behaviors of animals in a colony are relatively inflexible, and the relatively complex behaviors are controlled by the brain (Leadbeater and Chittka, 2007), we speculated that caste determination might be accompanied by differentiation of specific cell types. To identify the difference between queen and worker castes, we performed transcriptomic analysis of different caste brains at the single-cell level, and we identified a specific subtype of ensheathing cell with high levels of vg expression in queen brains. In particular, vg was more highly expressed compared to the average vg expression level in the ensheathing cells of workers. Importantly, we found that the knockdown of vg at the early larval stage, but not the late stage, significantly suppressed the developmental trajectory to queens. Therefore, our study not only identifies the high expression of vg in a subset of glial cells as a "molecular signature," but it also suggests that vg is involved in caste differentiation.
By analyzing transcriptomic features in this study, we obtained an important finding that glial cell populations vary significantly between queen and worker castes. Specifically, we found that the ratio of ensheathing glia was most variable among the castes and subcastes. The ensheathing glia population has been reported to be associated with caste transition, injury, and aging in ants (Sheng et al., 2020). In addition, the function of the ensheathing glia has been extensively studied in other model organisms, including Drosophila and mammals (Freeman, 2015;Ramon-Cueto and Avila, 1998). Our findings further emphasize the importance of this type of cell in regulating development and behavior.
Vitellogenin is a well-studied factor that is involved in a variety of biological processes, including reproduction, immunity, aging, and behavioral transition (Amdam et al., 2003;Havukainen et al., 2011;Ihle et al., 2015;Nelson et al., 2007;Seehuus et al., 2007). In particular, vg is synthesized by the fat body and then released to the hemolymph. Its expression in the hemolymph of queens has been reported to be much higher than that in workers (Barchuk et al., 2002;Guidugli et al., 2005). In this study, we identified a specific subtype of ensheathing glial cells that expresses vg at extremely high levels, suggesting that vg expression might be associated with the differentiation of this cell subtype. After knocking down vg with RNAi at an early larval stage, we found that the fate of queens was inhibited. In conclusion, our work identifies a subset of glial cells with vg high expression as a "biomarker" for the queen caste. It also suggests that vg plays a role in regulating caste differentiation in honeybees.

Limitations of the study
There are some limitations to this study. First, from a technical standpoint, the single-cell atlas cannot completely cover all the cell types, and rare cell populations might be omitted. Second, available genomic information of honeybees is insufficient and has not been fully elucidated. As a result, $30% of captured cells in our atlases are functionally undefined. Third, the variation between replicates increased during reclustering analysis, particularly for glia subtypes, which was also an inherent issue of the approach of singlecell RNA-seq. It is easy to produce the variation when we capture those rare cell types. Nevertheless, the findings from our ''FISH'' assays provided experimental evidence to support our conclusions.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  Figure S4 and

DECLARATION OF INTERESTS
The authors declare no competing interests.  (B-D) Simultaneous detection of tsf1 mRNA (magenta), idgf4 mRNA (green), vg mRNA (red), and (4 0 , 6-diamidino-2-phenylindole) DAPI-stained nuclei (blue), using confocal microscopy. The co-expression of tsf1 and idgf4 was used to mark ensheathing glia. Brain sections from the queen (B), nurse (C), and forager (D), respectively. Left: whole brain; Right: successively amplified regions with positive probe signals. Scale bars were indicated on the corresponding images. (E) Scheme of RNAi treatment assay based on the in-dish queen rearing system. First instar larvae were transferred into a 24-well plate and fed high nutritional foods. The double-strand RNAs were added or not added to the food in the treatment or control group, respectively. The pupae were transferred to another 24-well plate and kept still for eclosion. (F) The adults reared in the in-dish-rearing system are divided into three typical phenotypes according to birth weight. (G) Statistics for three phenotypes of four experimental groups (blank, ds GFP, ds vg and ds Tor). The phenotypes are colored according to different birth weight standards. The number of individuals and percentage of each phenotype are shown. See also Figure S5 and Table S12. ll OPEN ACCESS Leadbeater, E., and Chittka, L. (2007). Social learning in insects-from miniature brains to consensus building. Curr. Biol. 17, R703-R713. https://doi.org/10. 1016/j.cub.2007.06.012. Li, B., and Dewey, C.N. (2011 iScience Article iScience Article for 30 min. Washed for 5 min twice in PBST, and sections were dehydrated with gradient ascending methanol (25%, 50%, 75%, 100%) (5 min for each concentration on ice) when sections were treated in 100% meth-anol，-20 C overnight. Day2: Dehydrated sections were treated with gradient descending methanol (100%, 75%, 50%, 25%) to rehydrated, PBST washed for 5 min twice on ice. And proteinase K (10 mg/mL) for 3 min, washed with PBST 5 min twice, then fixed again with 4% paraformaldehyde for 1 h. Washed with PBST for 5 min thrice. And they were treated with wash buffer (wash buffer A: PBST = 1:1) for 10 min, wash buffer A (6 mL Stellarisâ RNA FISH Wash buffer A, 21 mL DEPC H2O, 3mL Formamide) for 5 min twice, HYB-(10 mL 20x saline sodium citrate, 20 mL formamide, 10 mL DEPC H2O, 40 mL Tweenâ20) for 5 min at room temperature then 1h at 37 C. Add probes (1: 600) diluted with HYB-(Probes for FISH were shown in Table S3). Placed sections at 37 C in an incubator for 18 h. Day3: Probes were removed thoroughly, HYB-30 min at 37 C, PBST 2 min twice, DAPI (0.5 mg/mL) 30 min at room temperature, PBST 10 min thrice. The sections were mounted with Antifade Mounting Medium. Image reconstruction was processed using confocal microscopy.

Confocal imaging
The whole-brain imaging was performed using Zeiss LSM 980 confocal microscopes with Plan-Apochromat 403/1.3 oil objective. Z-stack and tile scan features were used to image the large, wavy surfaces of the brain slices. The resulting tiles were then stitched into a single large image (ZEN 3.1 Blue software, Zeiss), which enabled visualization of the whole brain at high resolution. Imaris software (Imaris 9.0.1, Bitplane) was used to visualize images in 3D.
dsRNA synthesis and RNAi treatment dsRNAs were synthesized in an in vitro transcription reaction using T7 RiboMAXä Express Large-Scale RNA Production System (Promega). The DNA templates ($500bp) of dsRNAs were amplified from the queens' ovaries cDNA. The T7 promoter sequence was added to the 5 0 end of all the primers for DNA template amplification. The primers used in this study were listed in Table S1.
To promote the efficiency and specificity of gene silencing, three dsRNAs targeted to different gene regions were mixed in equal proportions and diluted to a final concentration of 1 mg/mL of total dsRNA. For one treatment on one 24-well plate filled with larvae, 2.4 mg dsRNA mixture and 8mL transfection polymer material (Poly[oxy-1,4-butanediyloxy(1-oxo-1,3-propanediyl)[[2-(dimethylamino)ethyl]imino](3-oxo-1,3-propanediyl)]) were diluted in 80 mL 25 mM NaOAc buffer (pH5.0), vortexed vigorously for 30 s and then incubated at room temperature for 20 min. Dropwise added 3 mL mixture to the region near the larva's head. The tips should avoid touching larva. This treatment was begun on the day after 1st instar larvae transfer and repeated once a day until the larval stage ended.

Synthesis of transfection polymer
To a solution of 1,4-butanediol diacrylate (10 mmol, 1982.2 mg, 1.982 mL) in THF (20 mL), diamine (10 mmol, 881.5 mg, 1.092mL) was added. The reaction mixture was stirred at 50 C for 48 h. After that, the reaction was cooled to room temperature and concentrated under reduced pressure. The polymer was prepared with DMSO into a solution of 500 mg/mL.

Quantitative real-time PCR
Total RNA was isolated with Trizol Reagent (Invitrogen). cDNA was synthesized using TransScriptâ First-Strand cDNA Synthesis SuperMix (TransGen). Quantitative RT-PCR was performed using AceQ Universal SYBR mixture (Vazyme) in duplicate on LightCycler 480 II (Roche). Template concentrations were normalized to endogenous reference RpL32. Primers for quantitative RT-PCR were designed by PerlPrimer and shown in Table S2.

10x data preprocessing
The 10x honeybee brain samples were each processed by the Cell Ranger (version 3.1.0) count pipelines provided by 10X Genomics to generate single feature-barcode matrices. The reference genome was built based on A. mellifera (assembly Amel_HAv3.1) from NCBI.

ll
OPEN ACCESS iScience 25, 104643, July 15, 2022 iScience Article Data processing by seurat The feature-barcode matrices of each sample (two biological replicates were performed for each queen, forager, and nurse experiment) were processed by Seurat (version 3.1.5) (Butler et al., 2018). Genes expressed in at least one cell were reserved for further processing. The duplicates of three types of honeybee brain samples were merged by the Seurat merge function. In addition, cells with unique feature counts between 200 and 5000, the number of UMIs between 500 and upper limit (upper limit = Q3 + 1.5 * (Q3 -Q1), Q1 and Q3 are the first and third quartiles) and less than 15% of mitochondrial genes were screened for further analysis.
After the filtered datasets were normalized using the default parameters, variable genes were found using the FindVariableFeatures function in Seurat with the parameters including selection.method = ''vst'', nfeatures = 3000. Then, three types of honeybee brain datasets were integrated using the IntegrateData function in Seurat with dims = 1:30 parameter.
Then, the data were reduced and clustered using principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) in Seurat. And Clustering of the data was processed using the FindClusters function with resolution = 1 parameter. Marks genes for each cluster were determined by Seurat FindAllMarkers function with parameters including only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.5, return.thresh = 0.01, test.use = ''bimod''.
For reclustering of glia and Kenyon cells, cells set were selected by their respective marker genes, and analysis was carried out by the above method. The only difference was that the datasets of glia and Kenyon cells adopt 0.2 and 0.5 resolutions respectively. The information of all described genes in this study is summarized in Table S13.

Analysis of bulk RNA-seq
Bulk RNA-seq libraries were paired-end-sequenced. FastQC software (version 0.11.9) was used for quality control. Adapters were trimmed off with the Trimmomatic program (version 0.39) (Bolger et al., 2014). After filtering, the paired reads were aligned to the honeybee reference genome (assembly Amel_HAv3.1) using STAR software (version 2.7.9a) (Dobin et al., 2013). The STAR mapping produced reads were accurately quantified by the RSEM program (version 1.3.1) (Li and Dewey, 2011) to obtain count matrices about gene expression. The count matrices were processed by DEseq2 (version 1.32.0) (Love et al., 2014) R package.

Analysis of neuropeptides
The expression profiles of 27 reported bee neuropeptides were obtained from single cell RNA-seq and bulk RNA-seq. TBtools software (version 1.082) (Chen et al., 2020) was used to generate a heatmap.

Differentially expressed genes
For single cell RNA-seq data, differentially expressed genes were calculated by Seurat FindMarkers function with parameters including logfc.threshold = 0, min.pct = 0, test.use = ''bimod''. The significance threshold was adjusted p value < 0.01, log2 foldchange of the average expression between the two groups >0.5 or < -0.5. For bulk RNA-seq data, adjusted p value < 0.01 and log2 foldchange >1 or < -1 in DEseq2 results were considered as differentially expressed genes.
Comparisons to bulk RNA-seq and single-cell RNA-seq Pseudo-bulk normalized expressions of single-cell RNA-seq were obtained from three Seurat objects of forager, nurse, and queen. Compared with the normalized counts of bulk RNA-seq, the genes that were simultaneously detected were retained. The normalized expressions of single-cell RNA-seq and bulk RNA-seq were converted by log2 to calculate Spearman's correlation coefficients using the stat_cor function from the ggpubr package.

Gene ontology enrichment analysis
The FindMarkers function of the Seurat R package was used to conduct data analysis. Genes with average log2 foldchange >0.5 and adjusted p value < 0.01 were considered as DE-Gs. GO term enrichment was performed using DAVID (Huang da et al., 2009a, b